DimPlot(obj.tmp, group.by = "orig.ident", pt.size = 0.2) + theme_bw() +
  xlim(c(11,13)) +ylim(c(-5,2)) +ggtitle("Cluster7:Dendritic Cells")

orig.ident_fils = paletteer::scale_fill_paletteer_d("ggsci::nrc_npg") 

res = "orig.ident"
obj.tmp@meta.data %>% ggplot(aes(!!sym(res), fill=!!sym(res))) + 
  geom_bar(alpha=0.7, color="grey5", size=0.1) +
  geom_text(stat="count", aes(label= ..count..), vjust=-0.5, size=3) +
  orig.ident_fils + 
  xlab("") + 
  theme_classic() +
  theme(legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust=1)) 

res = "orig.ident"
df = obj.srt@meta.data %>% dplyr::select(!!sym(res),RNA_snn_res.0.2)
df_all = df %>% dplyr::select(orig.ident) %>% table() %>% as.data.frame()  
df_7= df %>% filter(RNA_snn_res.0.2 == 7) %>% dplyr::select(orig.ident) %>% table() %>% as.data.frame()  

df_all$DC= df_7$Freq
df_all$DC_fraction= df_7$Freq/df_all$Freq *100
colnames(df_all)[2] = "Number_of_all_cells"

orig.ident_fils = paletteer::scale_fill_paletteer_d("ggsci::nrc_npg") 

ggplot(df_all, aes(x = orig.ident, y = DC_fraction, fill = orig.ident)) +
  geom_col(alpha = 0.7, color = "grey5", size = 0.1) +
  orig.ident_fils +
  geom_text(aes(label = round(DC_fraction, 2)), vjust = -0.5, size = 3) +  # 소수점 둘째자리까지 반올림하여 표시
  theme_classic() +
  theme(legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylab("Dendritic Cell Fraction in the sample") +
  xlab("")

# Define funtion 
#id1 = "treatment"
#id2 = "control"
#logfc = 0
deg.two.groups = function(obj.srt = obj.srt, treatment, control,logfc=0){
  Idents(obj.srt) = 'orig.ident'
  markers <-FindMarkers(
    obj.srt,
    logfc.threshold = logfc,
    ident.1 = treatment,
    ident.2 = control, slot= "data")
}

deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "P30342_D", 
                    control = "P30342_M")
deg1 =deg 
t= paste0(paste0("P30342_", "D"), "/",paste0("P30342_", "M"))

deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))

up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

deg %>% DT::datatable(extensions = "Buttons", options = list(dom="Bfrtip", buttons=c("csv","excel"), pageLength=10))

# GSEA related functions 
library(clusterProfiler)

perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    # Check the name of log2fc related 
    if ("avg_log2FC" %in% names(res)) {
      df <- res$avg_log2FC
    } else if ("log2FoldChange" %in% names(res)) {
      df <- res$log2FoldChange
    } else {
      stop("Neither avg_log2FC nor log2FoldChange columns found in the data frame.")
    }
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  result <- x@result %>% arrange(desc(NES))
  result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
  return(result)
}

# GSEA Plot 
gsea_nes_plot <- function(gsea.res, title, color="pvalue") {
  gsea.res = gsea.res %>% mutate(sig=ifelse(pvalue <= 0.05, "p value <= 0.05", "p value > 0.05"))
  # basic plot
  p <- gsea.res %>%
    ggplot(aes(reorder(ID, NES), NES)) +
    geom_col(aes(fill=!!sym(color)), color="grey5", size=0.15, alpha=0.8) +
    coord_flip() +
    labs(x="Pathway", y="Normalized Enrichment Score", title="GSEA") +
    theme_classic() +
    theme(axis.text.x = element_text(size=5, face = 'bold'),
          axis.text.y = element_text(size=6, face = 'bold'),
          axis.title = element_text(size=10)) +
    ggtitle(title)
  
  # color by color input type
  if (color == "pvalue") {
    p <- p + scale_fill_gradient(low = 'orangered', high = '#E5E7E9')
  } else if (color == "sig") {
    p <- p + scale_fill_manual(values = c("orangered", "#E5E7E9"))
  }
  return(p)
}

# Pathway selection : HALLMARK 
hallmark <- msigdbr::msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol) 

gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")

perform_GSEA2 <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$avg_log2FC
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  # Initialize .Random.seed object
  set.seed(NULL)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  
  return(x)
}

# Application 
gsea.res.raw = perform_GSEA2(res = deg1, ref = hallmark) 

deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "P30342_I", 
                    control = "P30342_M")
deg2 =deg 
t= paste0(paste0("P30342_", "I"), "/",paste0("P30342_", "M"))

deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))

up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

deg %>% DT::datatable(extensions = "Buttons", options = list(dom="Bfrtip", buttons=c("csv","excel"), pageLength=10))

gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")

perform_GSEA2 <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$avg_log2FC
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  # Initialize .Random.seed object
  set.seed(NULL)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  
  return(x)
}

# Application 
gsea.res.raw = perform_GSEA2(res = deg2, ref = hallmark) 

deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "P30342_A", 
                    control = "P30342_M")
deg3 =deg 
t= paste0(paste0("P30342_", "A"), "/",paste0("P30342_", "M"))

deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))

up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

deg %>% DT::datatable(extensions = "Buttons", options = list(dom="Bfrtip", buttons=c("csv","excel"), pageLength=10))

gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")

perform_GSEA2 <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$avg_log2FC
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  # Initialize .Random.seed object
  set.seed(NULL)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  
  return(x)
}

# Application 
gsea.res.raw = perform_GSEA2(res = deg3, ref = hallmark) 

deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "P30342_I", 
                    control = "P30342_D")
deg4 =deg 
t= paste0(paste0("P30342_", "I"), "/",paste0("P30342_", "D"))

deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))

up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

deg %>% DT::datatable(extensions = "Buttons", options = list(dom="Bfrtip", buttons=c("csv","excel"), pageLength=10))

gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")

perform_GSEA2 <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$avg_log2FC
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  # Initialize .Random.seed object
  set.seed(NULL)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  
  return(x)
}

# Application 
gsea.res.raw = perform_GSEA2(res = deg4, ref = hallmark) 

deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "P30342_A", 
                    control = "P30342_D")
deg5 =deg 
t= paste0(paste0("P30342_", "A"), "/",paste0("P30342_", "D"))

deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))

up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

deg %>% DT::datatable(extensions = "Buttons", options = list(dom="Bfrtip", buttons=c("csv","excel"), pageLength=10))

gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")

perform_GSEA2 <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$avg_log2FC
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  # Initialize .Random.seed object
  set.seed(NULL)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  
  return(x)
}

# Application 
gsea.res.raw = perform_GSEA2(res = deg5, ref = hallmark) 

deg = deg.two.groups(obj.srt = obj.tmp,
                    treatment = "P30342_A", 
                    control = "P30342_I")
deg6 =deg 
t= paste0(paste0("P30342_", "A"), "/",paste0("P30342_", "I"))

deg %>% ggplot(aes(avg_log2FC, -log10(p_val))) + 
  geom_point(size=0.1) +
  geom_vline(xintercept = 0, size=0.1) +
  theme_bw() + ggtitle(t)

deg = deg %>% mutate(DE= ifelse(p_val <= 0.05 & avg_log2FC >= log2(1.2), "UP",
                                ifelse(p_val <= 0.05 & avg_log2FC <= -log2(1.2), "DN",
                                "no significant")))
deg$DE = factor(deg$DE, levels = c("UP","DN","no significant"))

up = nrow(deg[deg$DE == "UP", ])
dn = nrow(deg[deg$DE == "DN", ])
deg %>% ggplot(aes(avg_log2FC, -log10(p_val), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(1.2), log2(1.2)), size=0.1, color="grey88") +
  geom_hline(yintercept = -log10(0.05), size=0.1) +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

deg %>% DT::datatable(extensions = "Buttons", options = list(dom="Bfrtip", buttons=c("csv","excel"), pageLength=10))

gsea.res = perform_GSEA(res = deg, ref = hallmark, pvalueCutoff = 1)

gsea_nes_plot(gsea.res = gsea.res, title = t, color = "pvalue") 

gsea_nes_plot <- function(gsea.res, title, color="pvalue") {
  gsea.res = gsea.res %>% mutate(sig=ifelse(pvalue <= 0.05, "p value <= 0.05", "p value > 0.05"))
  
  # 기본 플롯
  p <- gsea.res %>%
    ggplot(aes(reorder(ID, NES), NES)) +
    geom_col(aes(fill=!!sym(color)), color="grey5", size=0.15, alpha=0.8) +
    coord_flip() +
    labs(x="Pathway", y="Normalized Enrichment Score", title="GSEA") +
    theme_classic() +
    theme(axis.text.x = element_text(size=5, face = 'bold'),
          axis.text.y = element_text(size=6, face = 'bold'),
          axis.title = element_text(size=10)) +
    ggtitle(title)
  
  # color에 따른 색상 설정
  if (color == "pvalue") {
    p <- p + scale_fill_gradient(low = 'orangered', high = '#E5E7E9')
  } else if (color == "sig") {
    # sig 값이 "p value <= 0.05"가 없으면 기본 색상을 "#E5E7E9"로 설정
    if (!any(gsea.res$sig == "p value <= 0.05")) {
      p <- p + scale_fill_manual(values = c("#E5E7E9"))
    } else {
      p <- p + scale_fill_manual(values = c("orangered", "#E5E7E9"))
    }
  }
  return(p)
}
gsea_nes_plot(gsea.res = gsea.res, title = t, color = "sig")